options(warn=-1)
library(numDeriv)

###################
## AK1 Estimator
###################
AK1logLik <- function(para, mydata, z) {
  n = nrow(mydata);
  tauhat <- para[1];
  betap <- para[2];
  Coeff <- as.matrix(para[3:length(para)])
  if(ncol(mydata)>3){
    err <- mydata[,2] - as.matrix(mydata[,c(4:ncol(mydata))])%*%Coeff;
    esthat <- as.matrix(mydata[,c(4:ncol(mydata))])%*%Coeff;
  }else{
    err <- mydata[,2] - Coeff;
    esthat <- Coeff;
  }
  se <- mydata[,3];  
  t <- mydata[,2]/mydata[,3]
  cutoffs <- c(-1.96, 1.96)
  phat <- (abs(t)<1.96)*betap + (abs(t)>=1.96)*1
  meanbeta=rep(0,n);
  for (i in 1:n) {
    prob_mid <- (pnorm((cutoffs[2]*se[i]-esthat[i])/sqrt((tauhat)^2 +se[i]^2))
                 -pnorm((cutoffs[1]*se[i]-esthat[i])/sqrt((tauhat)^2 +se[i]^2)));
    prob_ex <- 1-prob_mid;
    meanbeta[i] <- betap*prob_mid + 1*prob_ex;
  }
  fX = dnorm(err, 0, sqrt(se^2 + tauhat^2));
  L <- (phat/meanbeta)*fX;
  logL <- log(L);
  LLH <- -sum(log(L));
  if(tauhat<0 | betap<0){ LLH <- 10^10; }
  if(z=="est"){return(LLH);}else{return(logL);}
}
AK1logLik_LLH <- function(para) {
  return(AK1logLik(para, AKdata, "est"))
}
ClusteredStdErr1 <-function(est, my.data, hessian) {
  g <- matrix(0, nrow=nrow(my.data), ncol=length(est))
  for(i in 1:length(est)){
    grad1 <- rep(0, length(est));
    grad2 <- rep(0, length(est));
    stepsize <- 10^(-6)
    grad1[i] <- -stepsize
    grad2[i] <- stepsize
    g[,i] <- (AK1logLik(est+grad2, my.data, "gra")-AK1logLik(est+grad1, my.data, "gra"))/(2*stepsize)
  }
  cluster_index<-my.data[,1]
  I<-order(cluster_index);
  cluster_index<-sort(cluster_index);
  g = g[I,]
  g = g - matrix(rep(apply(g,2,mean),length(I)),nrow=length(I),byrow=TRUE);
  gsum = apply(g,2,cumsum);
  index_diff <- cluster_index[-1] != cluster_index[-length(cluster_index)];
  index_diff <- c(index_diff,1);
  gsum = gsum[index_diff==1,]
  gsum=rbind(gsum[1,], diff(gsum));
  Sigma=1/(dim(g)[1]-1)*(t(gsum)%*%gsum);
  RobuClusteredStdErr<-diag(sqrt(nrow(my.data)*solve(hessian)%*%Sigma%*%solve(hessian)))
  return (RobuClusteredStdErr)
}
AK1Estimator <-function(AKdata) {
  InitialValue<-rep(1/2, (ncol(AKdata)-3))
  Result<-nlminb(AK1logLik_LLH, start=c(1,1,InitialValue),lower=c(0, 0, rep(-Inf,length(InitialValue))),upper=c(Inf, Inf, rep(Inf,length(InitialValue))), hessian=TRUE, control=list(iter.max=1000, abs.tol=10^(-20), eval.max=1000))
  EstCoefficients <- t(as.matrix(Result$par))
  colnames(EstCoefficients) <-  c("tau", "Betap", colnames(AKdata)[4:ncol(AKdata)])
  ClusStdErrors <- t(ClusteredStdErr1(Result$par, AKdata, hessian(AK1logLik_LLH, Result$par)))
  colnames(ClusStdErrors)<-  c("tau", "Betap", colnames(AKdata)[4:ncol(AKdata)])
  PValue <- dt(EstCoefficients/ClusStdErrors, df=nrow(AKdata)-length(EstCoefficients)-2)
  EstResults <- rbind(EstCoefficients, ClusStdErrors, PValue)
  rownames(EstResults) <- c("Coefficients", "std.Err.Clustered", "p-value")
  AKEstimationResults<-list(Est.Type = 'AK1 (symmetric)',
                            nlmOutput = Result,
                            EstResults = EstResults,
                            LogLiklihood = -Result$objective,
                            SampleSize = nrow(AKdata),
                            aic =  2*ncol(EstResults) - 2*(-Result$objective), 
                            bic =  log(nrow(AKdata))*ncol(EstResults) - 2*(-Result$objective))
  return(AKEstimationResults)
}
###################


###################
## AK2 Estimator
###################
AK2logLik <- function(para, mydata, z) {
  n = nrow(mydata);
  tauhat <- para[1];
  beta1 <- para[2];
  beta2 <- para[3];
  beta3 <- para[4];
  Coeff <- as.matrix(para[5:length(para)])
  if(ncol(mydata)>3){
    err <- mydata[,2] - as.matrix(mydata[,c(4:ncol(mydata))])%*%Coeff;
    esthat <- as.matrix(mydata[,c(4:ncol(mydata))])%*%Coeff;
  }else{
    err <- mydata[,2] - Coeff;
    esthat <- Coeff;
  }
  se <- mydata[,3];  
  t <- mydata[,2]/mydata[,3]
  cutoffs <- c(-1.96, 0, 1.96);
  phat <- (t<=-1.96)*beta1 + (t>-1.96 & t<=0)*beta2 + (0< t &t<1.96)*beta3 + (t>=1.96)*1
  meanbeta=rep(0,n);
  for (i in 1:n) {
    prob_vlow <- pnorm((cutoffs[1]*se[i]-esthat[i])/sqrt((tauhat)^2 +se[i]^2))
    prob_low <- (pnorm((cutoffs[2]*se[i]-esthat[i])/sqrt((tauhat)^2 +se[i]^2))
                 -pnorm((cutoffs[1]*se[i]-esthat[i])/sqrt((tauhat)^2 +se[i]^2)));
    prob_upper<- (pnorm((cutoffs[3]*se[i]-esthat[i])/sqrt((tauhat)^2 +se[i]^2))
                  -pnorm((cutoffs[2]*se[i]-esthat[i])/sqrt((tauhat)^2 +se[i]^2)));
    prob_vupper<- 1-prob_vlow-prob_low-prob_upper;
    meanbeta[i] <- beta1*prob_vlow + beta2*prob_low + beta3*prob_upper + 1*prob_vupper;
  }
  fX = dnorm(err, 0, sqrt(se^2 + tauhat^2));
  L <- (phat/meanbeta)*fX;
  logL <- log(L);
  LLH <- -sum(log(L));
  if(tauhat<0 | beta1<0 | beta2<0 | beta3<0 ){ LLH <- 10^10; }
  if(z=="est"){return(LLH);}else{return(logL);}
}
AK2logLik_LLH <- function(para) {
  return(AK2logLik(para, AKdata, "est"))
}
ClusteredStdErr2 <-function(est, my.data, hessian) {
  g <- matrix(0, nrow=nrow(my.data), ncol=length(est))
  for(i in 1:length(est)){
    grad1 <- rep(0, length(est));
    grad2 <- rep(0, length(est));
    stepsize <- 10^(-6)
    grad1[i] <- -stepsize
    grad2[i] <- stepsize
    g[,i] <- (AK2logLik(est+grad2, my.data, "gra")-AK2logLik(est+grad1, my.data, "gra"))/(2*stepsize)
  }
  cluster_index<-my.data[,1]
  I<-order(cluster_index);
  cluster_index<-sort(cluster_index);
  g = g[I,]
  g = g - matrix(rep(apply(g,2,mean),length(I)),nrow=length(I),byrow=TRUE);
  gsum = apply(g,2,cumsum);
  index_diff <- cluster_index[-1] != cluster_index[-length(cluster_index)];
  index_diff <- c(index_diff,1);
  gsum = gsum[index_diff==1,]
  gsum=rbind(gsum[1,], diff(gsum));
  Sigma=1/(dim(g)[1]-1)*(t(gsum)%*%gsum);
  RobuClusteredStdErr<-diag(sqrt(nrow(my.data)*solve(hessian)%*%Sigma%*%solve(hessian)))
  return (RobuClusteredStdErr)
}
AK2Estimator <-function(AKdata) {
  InitialValue<-rep(1/2, (ncol(AKdata)-3))
  Result<-nlminb(AK2logLik_LLH, start=c(1,1,1,1,InitialValue), lower=c(0, 0, 0, 0, rep(-Inf,length(InitialValue))),upper=c(Inf, Inf, Inf, Inf, rep(Inf,length(InitialValue))), hessian=TRUE, control=list(iter.max=1000, abs.tol=10^(-20), eval.max=1000))
  EstCoefficients <- t(as.matrix(Result$par))
  colnames(EstCoefficients) <-  c("tau", "Beta1", "Beta2", "Beta3", colnames(AKdata)[4:ncol(AKdata)])
  ClusStdErrors <- t(ClusteredStdErr2(Result$par, AKdata, hessian(AK2logLik_LLH, Result$par)))
  colnames(ClusStdErrors)<-  c("tau", "Beta1", "Beta2", "Beta3", colnames(AKdata)[4:ncol(AKdata)])
  PValue <- dt(EstCoefficients/ClusStdErrors, df=nrow(AKdata)-length(EstCoefficients)-2)
  EstResults <- rbind(EstCoefficients, ClusStdErrors, PValue)
  rownames(EstResults) <- c("Coefficients", "std.Err.Clustered", "p-value")
  AKEstimationResults<-list(Est.Type = 'AK1 (asymmetric)',
                            nlmOutput = Result,
                            EstResults = EstResults,
                            LogLiklihood = -Result$objective,
                            SampleSize = nrow(AKdata),
                            aic =  2*ncol(EstResults) - 2*(-Result$objective), 
                            bic =  log(nrow(AKdata))*ncol(EstResults) - 2*(-Result$objective))
  return(AKEstimationResults)
}
###################